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1  Introduction  and  Overview 


While  computational  methods  have  become  increasingly  refined  and  accurate,  their  reliance  on  exact 
data,  e.g.,  complete  descriptions  of  geometries,  materials,  sources,  etc.,  is  emerging  as  a  bottleneck 
in  the  modeling  of  problems  of  realistic  complexity.  For  instance,  if  one  attempts  to  model  an 
experiment,  a  classic  computational  approach  would  require  knowledge  of  a  degree  of  detail  which 
is  unrealistic  and  often  impossible  to  obtain,  e.g.,  one  cannot  hope  to  control  all  elements  of  an 
experiment,  measure  all  details  of  an  initial  condition  or  geometry,  know  the  microstructure  of  all 
materials,  etc. 

The  usual  approach  to  dealing  with  this  lack  of  knowledge  or  uncertainty  is  to  simply  assume 
some  mean  parameters  and  compute  the  corresponding  solution.  If  the  solution  is  robust  to  pa¬ 
rameter  variation,  this  is  indeed  a  reasonable  approach.  However,  for  general  problems,  where  the 
sensitivity  of  parts  of  the  solution  can  be  significant,  a  solution  based  on  mean  parameters  is  not 
likely  to  match  very  well  with  experiments  and,  thus,  is  not  a  good  predictive  tool.  An  experimen¬ 
talist  would  naturally  deal  with  this  exact  problem  by  repeating  the  experiments  and,  subsequently, 
compute  not  only  mean  results  but  also  error  bars — reflecting,  at  least  partly,  the  sensitivity  of  the 
solution. 

It  is  natural  to  strive  to  achieve  similar  abilities  in  computational  modeling  effort,  e.g.,  compute 
solutions  or  measures  of  interest  and  their  associated  sensitivities.  Clearly,  if  a  particular  measured 
entity  is  highly  sensitive  there  is  no  reason  to  expect  that  this  matches  experiments  exactly.  Thus, 
we  would  like  to  be  able  to  model  the  impact  of  the  uncertainty,  assumed  to  have  certain  properties 
derived  from  experiments  or  otherwise,  on  the  computed  results,  essentially  resulting  in  an  ensample 
of  possible  solution  values  with  an  associated  probability.  This  would  immediately  enable  the 
computation  of  statistical  moments,  e.g.,  means  and  variances,  and  of  other  valuable  information 
about  the  sensitivity  of  solutions  and  derived  quantities. 

In  this  project  we  have  pursued  this  goal  and  we  present  here  a  systematic,  accurate,  and  efficient 
way  of  addressing  this  type  of  problem,  essentially  enabling  one  to  compute  with  an  ensample  of 
data  and,  subsequently,  obtain  a  full  space-time  ensample  of  solutions  with  an  associated  probability 
density.  It  is  important  to  realize  that  is  this  not  a  simple  situation,  since  solutions  may  vary 
nonlinearly  in  the  uncertainty  due  to  stochastic  correlations  even  if  the  deterministic  problem  is 
linear,  e.g.,  Maxwell’s  equations. 

A  standard  way  of  addressing  problems  of  the  type  mentioned  in  the  above  is  through  Monte 
Carlo  sampling,  e.g.,  run  a  deterministic  code  a  large  number  of  times  and,  subsequently,  extract 
the  statistics  of  interest.  The  main  problem  with  this  approach  is  the  very  slow  convergence  rate, 
(9(7V-1/2),  with  N  being  the  number  of  samples,  which  makes  even  the  computation  of  mean 
solutions  expensive  and  accurate  recovery  of  higher  moments,  e.g.,  variances,  prohibitive. 

Our  platform  on  which  to  demonstrate  this  approach  is  the  time-domain  Maxwell’s  equations, 
solved  using  a  high-order  accurate  discontinuous  Galerkin  method.  However,  the  basic  elements  of 
the  formulation  are  general  and  can  be  used  with  any  computational  kernel. 

We  have  focused  on  the  development  of  efficient  and  accurate  means  by  which  to  provide  quan¬ 
titative  measures  on  the  impact  on  the  radar  cross  section  by  uncertainties  in  general  scattering 
problems.  The  uncertainties  can  be  present  in  the  sources,  the  materials,  or  the  geometries. 

A  summary  of  the  main  accomplishments  are 

•  Validated  the  use  of  Wiener  chaos  expansion  methods  for  a  variety  of  Maxwell  benchmark 
problems,  including  internal,  external,  time-domain,  and  harmonic  problems. 
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•  Developed  a  collocation  form  of  the  chaos  expansion  method  to  dramatically  decrease  the 
computational  cost  of  using  such  methods. 

•  Developed  the  chaos  techniques  to  enable  the  modeling  to  uncertain  sources,  materials,  and 
geometries.  The  latter  has  resulted  in  the  introduction  of  a  new  stochastic  mapping  technique 
to  enable  the  modeling  of  general  scatterers. 

•  Developed  techniques  to  model  spatially  correlated  uncertainty  based  on  statistical  knowledge 
of  the  uncertainty. 

•  Developed  a  theoretical  understanding  of  the  convergence  rate  and  long  time  behavior  of  the 
chaos  expansion  methods. 

•  Demonstrated  the  ability  to  apply  such  techniques  to  quantify  the  impact  of  uncertainty  on 
scattering  problems  and  the  resulting  sensitivity  of  the  RCS. 

In  the  following  we  shall  elaborate  a  little  more  on  some  of  these  many  aspects. 

2  Maxwell’s  equations 

Let  us  consider  a  general  domain  D  and  let  E5  and  H5  denote  the  scattered  electric  and  magnetic 
fields,  respectively.  With  e(x)  and  //(x)  being  the  local  permittivity  and  permeability,  and  er(x)  the 
conductivity  of  the  media,  the  time-dependent  Maxwell’s  equations  in  the  scattered  field  formulation 
are  given  as 

n  -ps 

e-7—  =  V  x  Hs  +  cEs  +  SE,  (1) 

ot 

<9H,S 

/i-^-  =  —  V  xE’  +  SH ,  (2) 

where,  as  is  common  in  time-domain  schemes,  we  have  neglected  divergence  constraints,  assuming 
that  these  amount  to  constraints  on  the  initial  conditions. 

The  source  terms,  Sh  and  S^,  appearing  on  the  right-hand  side  of  ( 1)— (2) ,  take  the  form 


SB  =  -(e-e^  +  ^-^E' 

(3) 

s  =  </*/*>*. 

(4) 

where  the  incident  field  (E*,Hl)  is  a  solution  to  Maxwell’s  equations  in  a  media  of  permittivity, 
permeability,  and  conductivity — £l(x),  /i*(x),  crl(x),  respectively. 

We  now  turn  our  attention  to  boundary  conditions.  Along  a  perfect  electric  conductor  (PEC), 
tile  boundary  conditions  on  the  electric  field  are 

n  x  E'  =  0, 

(5) 

where  n  is  the  outward  pointing  normal  vector  at  the  surface  and  E(  —  E' 
For  the  total  magnetic  field  H  ,  the  condition  is 

+  E 5  is  the  total  field. 

H'  n  =  0. 

(6) 
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Along  the  interface  of  two  dielectric  bodies,  denoted  by  the  subscripts  1  and  2,  we  have 

n  x  (EJ  -  E|)  =  0  and  n  x  (HJ  -  Us2)  =  0;  (7) 

i.e.,  all  tangential  components  are  continuous. 


3  Numerical  scheme  for  the  deterministic  case 


Before  turning  our  attention  to  the  case  including  uncertainties,  let  us  briefly  describe  the  computa¬ 
tional  methods  used  for  solving  Maxwell’s  equations  in  the  physical  space.  We  use  a  discontinuous 
Galerkin  formulation;  the  solution  will  be  discontinuous  between  elements.  This  offers  a  num¬ 
ber  of  advantages  over  widely  used  alternative  techniques,  e.g.,  geometric  flexibility  through  fully 
unstructured  grids,  high-order  accuracy  to  enable  accurate  wave  propagation  over  long  distance 
with  a  coarse  resolution,  and  very  high  computational  efficiency  through  explicit  time  stepping  and 
high  parallel  performance.  The  method  has  been  discussed,  analyzed,  and  validated  extensively  for 
Maxwell’s  equations  and  we  shall  simply  sketch  its  main  components. 

We  rewrite  Maxwell’s  equations  (l)-(2)  in  conservation  form 

+  V  •  F(q)  =  S,  (8) 

where  q  is  the  state  vector  formed  by  the  scattered  electric  field  and  the  magnetic  field 


q  = 


0) 


and  the  components  of  the  tensor  F  are  defined  by 


F,(q) 


-e2  xH\ 
e*  x  E  )' 


(10) 


where  e*  denotes  the  Cartesian  unit  vectors.  On  the  right-hand  side  of  (8),  S  =  [S^,  S^]  is  the 
source  term,  which  depends  on  the  incident  field,  and  the  material  matrix  Q  is  a  diagonal  matrix 
with  values  (e,  e,  e,  fi.  /i,  /i)  on  its  diagonal. 

We  assume  that  the  computational  domain,  ft,  is  tessellated  by  triangles  or  tetrahedrons,  D, 
and  we  represent  the  local  solution  q yv  as 


TV 

q/v(x,<)  =  ^qi(t)Li(x),  (11) 

1=1 

where  qz  reflects  nodal  values,  defined  on  the  element,  and  Li(x)  signifies  an  nth  order  Lagrange 
polynomial,  associated  with  grid  points. 

The  discrete  solution,  qv,  of  Maxwell’s  equations  is  required  to  satisfy 

In  +  ^  F(qjV)  ~  S/v)  L'(x)dx  =  '  [F(q/v)  -  F’jL^xjdx.  (12) 

In  (12),  F*  denotes  a  numerical  flux  and  n  is  an  outward  pointing  unit  vector  defined  at  the  boundary 
dD  of  the  clement  D.  Note  that  this  is  an  entirely  local  formulation,  and  relaxing  the  continuity 
of  the  elements  decouples  the  elements,  resulting  in  a  block-diagonal  global  mass  matrix  which  can 
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be  trivially  inverted  in  preprocessing.  The  price  to  pay  for  this  is  the  additional  degrees  of  freedom 
needed  to  support  the  local  basis  functions.  For  high-order  elements,  this  is,  however,  only  a  small 
fraction  of  the  total  number  of  degrees  of  freedom.  Given  the  linearity  of  Maxwell’s  equations,  and 
for  stability  reasons,  it  is  natural  to  use  an  upwinding  flux  which  takes  the  form 


n  •  [F(qjv)  -  F*] 


Z  'nx  (nx  [EN]-Z+[Hyv])  \ 
r'iix  (nx  [E/vJ-Z+fH/v])  /’ 


(13) 


where  the  bracket  [q]  =  q“  —  q+  denotes  the  jump  in  the  field  values  (q“  is  the  local  value  and 
qf  is  the  neighboring  value)  across  an  interface,  Z*  is  the  local  impedance,  and  Y±  is  the  local 
conductance, 


Z± 


1 


(14) 


Z  and  Y  in  (13)  are  the  sums 


z  =  z+  +  z~,  y  =  + 


(15) 


After  discretization  of  the  operators  and  evaluation  of  the  integrals  appearing  in  (12),  the  problem 
can  be  rewritten  in  matrix-vector  form 


QM-^-  4-  5  •  FjV  -  M SN  =  Fh  •  [F*  -  F*].  (16) 

The  matrices  A/,  5,  and  F  represent  the  local  mass-,  stiffness-,  and  face-integration  matrices. 
Note  that  the  local  nature  of  the  scheme  allows  for  the  use  of  an  explicit  solver  for  the  time 
discretization  of  (16).  This  is  traditionally  done  using  an  explicit  Rungc-Kutta  method,  although 
suitable  alternatives  arc  plentiful. 


4  Accounting  for  the  uncertainty 


To  deal  with  the  actual  lack  of  detailed  knowledge  leading  to  the  uncertainty,  we  must  make  some 
educated  guesses — a  model — about  the  nature  of  the  possible  variations  in  the  data.  These  models 
can  be  based  on  pure  speculation,  on  measured  data,  or  on  other  available  information.  Once  this 
is  done,  however,  we  must  introduce  this  into  the  computational  approach  in  an  efficient,  accurate, 
and  robust  manner. 

Consider,  as  a  simple  example,  the  wave  equation 


du  ,n,du 
-  +  a(0)— =  0, 


(17) 


where  0  is  a  random  parameter  with  some  associated  probability  density  function  (PDF).  Thus, 
this  represents  an  ensemble  of  wave  problems,  each  with  an  individual  phase  speed  occurring  with 
a  probability  given  by  the  PDF;  i.e.,  the  solution  u  is  not  only  a  function  of  (x,  t)  but  also  of  0 ,  and 
the  equation  is  a  stochastic  wave  equation. 

The  actual  form  of  a{6)  may  not  be  known,  causing  the  introduction  of  the  uncertainty.  However, 
it  seems  quite  reasonable  that  in  many  situations  one  can  measure  or  conjecture  the  average  of  a 
and  perhaps  its  variation;  i.e.,  one  can  reasonably  assume  that  a(0)  varies  in  a  certain  way  related 
to  a  normal  distribution  with  a  given  mean  and  variance. 
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The  question  remains,  however,  how  this  uncertainty  will  affect  the  solution,  u(x,£,  (9),  and, 
more  often,  its  moments  such  as  the  mean  and  the  variance.  This  is  not  a  trivial  question,  as  the 
uncertainty  essentially  renders  the  simple  linear  problem  considerably  more  complex  due  to  possible 
stochastic  dependence  between  a  and  u. 

One  could  naturally  solve  the  above  problem  for  a  large  number  of  values  of  9 ,  taken  from  a 
proper  distribution,  and,  subsequently,  compute  the  appropriate  moments.  This  is  the  essence  of 
a  Monte  Carlo  method  and  suffers  from  a  very  slowly  converging  result  as  1  /y/N,  with  N  being 
the  number  of  samples.  If  higher-order  moments,  e.g.,  the  variance  or  sensitivity  of  the  result,  are 
required,  this  is  often  prohibitive. 

In  the  above  example,  the  uncertainty  enters  through  the  phase  speed.  However,  as  we  shall  sec 
shortly,  dealing  with  uncertainty  in  initial  conditions,  boundary  conditions,  sources,  or  computa¬ 
tional  domain/geometries  can  be  done  in  an  entirely  similar  fashion. 

In  this  section,  we  shall  discuss  two  variations,  referred  to  as  stochastic  Galerkin  and  collocation, 
respectively,  of  the  same  underlying  result.  This  enables  one  to  discretize  stochastic  partial  differ¬ 
ential  equations  (SPDEs)  to  recover  systems  of  deterministic  problems  which  we  can  subsequently 
solve  with  the  methods  discussed  later.  As  we  shall  experience  through  examples,  both  techniques 
offer  excellent  accuracy  at  a  severely  reduced  computational  cost  as  compared  to  straightforward 
Monte  Carlo  methods. 


4.1  The  homogeneous  chaos  expansion 

The  key  result  on  which  we  shall  rely  is  due  to  Wiener  and  shows  that  the  chaos  expansion  can  be 
used  to  approximate  any  functional  in  L2(fl,P),  where  V  is  a  Gaussian  measure  on  ft.  For  such  a 
random  variable  X(9),  the  chaos  expansion  is  written  as 

X(9)  =  ci,qHo  (18) 

d 

+  £>,#1  (6,(0)) 

*1  =  1 

*1=1  *2  =  1 
d  ii  *2 

+  £  £  £  aili2i3  #3(6, (*).&(*), &3(*)) 
tl  =  l  12  =  1  *3=1 

+  •  *  •  , 


where  f  =  (£i(0),. . .  ,£<*(#))  represents  d  independent  Gaussian  variables  with  zero  mean  and  unit 
variance,  each  depending  on  the  random  event  9,  and  Hn  arc  the  multivariate  Hermite  polynomials 
defined  as 


e*«T<(— 1)" 


dn 


-if1 

D  O'4 


(19) 


The  number  of  terms  in  the  expansion  (18)  grows  as 


(n  +  d)! 
n\d\ 


(20) 


where  n  is  the  length  of  the  Hermite  expansion  and  d  is  the  dimension  of  the  Gaussian  random 
space.  The  polynomial  homogeneous  chaos  expansion  forms  a  complete  orthogonal  basis  in  the 
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space  of  Gaussian  variables,  i.e., 


(mn  (9),  •  •  •,&.(*)),  Mn  (0), . .  .,&,(*))>  (21) 

.  f  •  •  •  ’bn m e-itTtjr  _ 

]ra 

where  i  denotes  the  multi-index  (h, . . . ,  in)  and  i\  =  i\ in\.  Thus  there  is  an  intimate  relation 
between  the  Hermite  polynomials,  orthogonal  under  the  Gaussian  weight,  and  the  representation 
of  random  variables  taken  from  a  Gaussian  distribution.  One  way  of  interpreting  the  homogeneous 
chaos  expansion  is  that  a  general  random  variable  can  be  expressed  in  terms  of  simpler  Gaussian 
variables  for  which  we  can  construct  an  efficient  computational  approach.  Clearly,  if  the  random 
variable  is  far  from  Gaussian,  many  terms  in  the  expansion  will  be  needed;  i.e.,  n  must  be  large.  An 
alternative  is  to  use  an  expansion  in  terms  of  other  random  variables  with  an  associated  distribution 
closer  to  what  is  expected  at  input  or  output.  For  notational  convenience,  we  can  rewrite  (18)  in 
the  form 

p 

*(*)  =  £  WO,  (22) 

3- 1 

where  there  is  a  one-to-one  correspondence  between  the  functions  Hi  and  tyj  and  also  between  the 
coefficients  ailv..tjp  and  bj.  In  the  case  of  a  general  stochastic  process,  these  coefficients  will  be  t  ime 
dependent. 

To  model  the  impact  of  uncertainty  on  the  propagation  of  electromagnetic  waves,  we  include  the 
randomness  in  the  usual  spatial-temporal  dimensions;  i.e.,  the  electric  field  and  the  magnetic  field 
become  E(x,£,0)  and  H(x,£,0),  reflecting  that  the  fields  are  functions  of  d  independent  random 
variables,  (&,(0), . . .  ,&d(0))- 

In  the  following  we  shall  discuss  in  some  detail  how  this  can  be  utilized  to  construct  an  efficient 
computational  method.  For  simplicity  of  the  discussion,  we  assume  in  what  follows  that  one  Gaus¬ 
sian  variable  suffices  to  represent  the  process  (i.e.,  d  =  1).  However,  the  formulation  is  general  and 
applies  to  problems  requiring  many  random  variables  to  describe  the  stochastic  processes. 


4.2  Stochastic  Galerkin  formulation 

Using  the  chaos  expansion,  we  can  express  q(x,  t,  0)  =  (E(x,  t ,  6),  H(x,  t ,  0))T  as 

p 

q(x>  t,  9)  =  q'(xi  (23) 

i=  1 

We  can  write  the  computational  scheme,  which  takes  into  account  randomness  in  a  general  setting, 
as 

f  Q(0)M^  +  S  ■  Fn  -  MS(0)n  =  Fn  ■  [FN  -  F*], 
l  q/v(x,t  =  0,0)  =  f(x,0), 

where  the  initial  conditions  are  given  by  the  function  f  =  f(x.$)  and  the  unknown  vector  q,v  is 
given  by  (23).  As  a  first  step,  we  discretize  (24)  in  the  random  space  using  a  Galerkin  approach. 
Multiplying  (24)  by  a  test  function  4^(0),  replacing  q/v  by  its  chaos  expansion,  and  using  the  scalar 
product  defined  by  (21),  we  obtain 

p  do1  p 

vfc  G  [1,P]  :  ^(Q9i,9k)M-te  +  k\S-FkN-MSkN  =  F'52n-\EiN-Fi*],  (25) 

1=1  i=l 
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where  we  have  used  property  (21).  The  detailed  expression  of  the  terms  on  the  right-hand  side  is 


/  (z  n  x  n  x  [E*]  -  (Z+Z  *fc)  n  x  [IV]  \ 

\  (F'^^nxnx  +  [E*]  ) 


(26) 


Recall  that  Z  and  Y  may  depend  on  the  material  properties  and,  thus,  may  include  uncertainty/randomness. 

The  initial  conditions  in  (24)  also  need  to  be  projected  onto  the  chaos  basis  to  give  an  initial 
condition  for  each  mode  of  qlN  in  the  chaos  expansion,  i.e., 

•  Vie  [1,P]  :  q*v(x,t  =  0)  =  l(f(x,5),^).  (27) 

1 1 


Once  the  vectors  {qyv)i<i<P  of  the  system  (31)  have  been  computed,  we  have  available  at  every 
point  in  space  an  approximation  to  the  probability  density  of  the  solution  to  the  system. 

Considering  (31),  we  observe  that  we  have  managed  to  recast  the  general  stochastic  problem  into 
a  system  of  P  coupled  deterministic  problems  which  we  can  now  discretize  in  space/timc  exactly  as 
discussed  in  section  3 — or  in  any  other  preferred  way. 


4.3  Stochastic  collocation  formulation 

The  idea  of  the  stochastic  collocation  formulation  is  to  replace  the  expression  of  the  electric  and 
magnetic  fields  (23)  in  the  polynomial  chaos  expansion  with  a  Lagrangian  polynomial  basis;  i.e.,  we 
would  have 

p 

q(x,  t,e)  =  J2  q*(x,  t)Li(O),  (28) 

i=l 

where  forms  a  Lagrangian  polynomial  basis  of  degree  ( P  -  1).  In  the  deterministic 

community,  (23)  would  be  referred  to  as  a  modal  expansion  and  (28)  as  a  nodal  expansion.  Since 
6  in  (28)  is  a  Gaussian  variable  with  zero  mean  and  unit  variance,  it  seems  natural  to  use  Gauss- 
Hermite  collocation  points  as  the  basis  for  the  Lagrange  polynomials,  as  that  also  enables  accurate 
evaluation  of  integrals. 

The  stochastic  collocation  method  is  essentially  a  deterministic  or  lattice  Monte  Carlo  method 
with  samples  {\/20j}i<j<p.  The  crucial  difference  is  that  a  typical  Monte  Carlo  simulation  consist¬ 
ing  of  P  realizations  will  converge  at  the  slow  rate  1  /\[P,  whereas  the  convergence  of  the  stochastic 
collocation  method  is  much  faster,  as  will  be  shown  by  numerical  examples.  A  simple  analogy  is 
that  the  approximation  of  a  smooth  function,  the  PDF,  is  done  most  efficiently  by  representing 
it  by  smooth  polynomials  based  on  points  well  suited  for  interpolation.  In  a  simple  Monte  Carlo 
approach,  the  interpolation  points  are  random,  leading  to  a  poor  convergence  rate. 

We  presented  the  stochastic  collocation  method  for  a  random  space  of  size  d  =  1.  The  gener¬ 
alization  to  higher  random  spaces  is  straightforward,  using  tensor  products  of  quadrature  points. 
Note  that  for  higher  random  spaces  and  for  purposes  of  efficiency,  it  might  be  necessary  to  use 
sparse  grid  methods 


4.4  Modelling  of  a  random  surface 


We  consider  an  object  whose  shape  can  vary  in  a  random  fashion.  As  an  example,  Figure  1  shows  a 
disc  which  is  modelled  by  (n  -  1)  line  segments  [X*Xi+i]i<i<n_i  from  a  finite  element  mesh.  Those 
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Figure  1:  Points  X{  =  (r;  cos(0z),  rz  sin(0z))  defining  the  boundary  of  an  object 


line  segments  are  defined  by  the  points  of  polar  coordinates  Xz  =  (rz  cos (0*),  rz  sin(0z))i<j<n.  We  now 
assume  that  the  point  X7  can  be  moved  randomly  by  a  quantity  JX7;  =  ( Sri  cos (0*),  Jr*  sin(0j))i<j<n 
to  take  a  new  position  (Xz  -f  JXZ)  =  ((rz  +  Jrz)  cos(0*),  (rz  +  Jrz)  sin(0i))i<*<n.  We  further  assume 
that  two  points  (say  Xz  and  Xj)  with  polar  angles  Oi  and  close  to  each  other  should  have  a 
random  height  ||JX;||  close  to  ||JXj||.  This  is  done  by  introducing  the  covariance  matrix  K  such 
that: 


where 


f  Oi  if  9i  e  [0, 7r] 

\  7 t  -di  if  6i  e]0, 27r] 


(29) 


(30) 


In  the  definition  of  Pi  we  have  to  separate  the  cases  E  [0, 7r]  and  Oi  e]0,  27t]  to  ensure  that  two 
points  on  the  random  surface,  one  with  a  polar  angle  close  to  0  and  the  other  one  with  a  polar  angle 
close  to  27 r,  be  strongly  correlated.  In  the  relation  (29),  b  is  a  parameter  which  can  control  how 
correlated  two  points  X*  and  Xj  can  be  and  c  is  another  parameter  which  controls  the  roughness 
of  the  surface  (the  magnitude  of  ||£XZ||  for  i  =  1,  ..,n  is  directly  proportional  to  c).  Furthermore,  c 
is  the  standard  deviation  of  each  component  of  Sr. 


To  generate  a  random  surface,  the  problem  can  be  formulated  as  follows:  find  a  random  vector 
Sr  =  ( Sr\,Sr2 ,  ...,5rn)r  with  a  given  covariance  matrix  K  that  will  generate  the  new  random  surface 
(X;  -f  SXp  =  ((rz  +  Sri)  cos(Oi),  (ri  +  Sri)  sin(0i))i<i<n.  This  procedure  is  illustrated  in  Figure  2, 
where  a  triangle  (j4,Xj,Xj+i)  having  one  side  [X;Xz+i]  sitting  on  the  boundary  of  the  object  is 
distorted  into  a  new  triangle  (^4,  Xz  +  <5Xj,  Xi+i  +  £Xz+i)  which  side  [(Xz  +  <5Xz)(Xz+i  +  <5X;+i)]  will 
compose  the  new  shape  of  the  object.  From  a  practical  point  of  view,  this  can  be  easily  done  by  first 
generating  a  vector  G  =  \/3(ai,  Q2,  ctn)T ,  where  oi,a2j  •••>  are  random  independent  variables 
following  (for  example)  a  uniform  law  on  the  interval  [-1,1].  Since  K  is  symmetric  definite  positive, 
it  can  be  decomposed  as  K  =  PDP"1,  where  D  is  a  diagonal  matrix  with  positive  eigenvalues  on  its 
diagonal.  Then,  it  can  be  shown  that  the  vector  Jr  given  by  Jr  =  PD1/2G  will  be  a  random  vector 
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A 


Figure  2:  Illustration  showing  how  the  points  on  the  boundary  of  an  object  are  moved  randomly. 

with  covariance  matrix  K.  In  terms  of  implementation,  one  generates  independents  pseudo-random 
numbers  a i,a2,  ...,an  with  a  uniform  law  on  the  interval  [—1, 1],  then  the  vector  Jr  is  given  by  the 
matrix- vector  product  PD1  2G  (the  matrix  PD1//2  can  be  pre-computed  and  stored  once  for  all). 
The  coordinates  of  the  points  of  the  new  random  surface  will  be  (X;  +  SXt)  for  2=1,  ...,n.  The 
procedure  described  here  for  a  disc  can  be  easily  adapted  to  objects  with  more  general  shapes,  as 
will  be  shown  in  the  numerical  experiments  section. 

As  an  example,  we  have  considered  a  somewhat  simplified  rocket  in  Figure  3.  Figure  4  shows 
a  zoom  of  its  front  part  for  the  original  (non-distorted)  rocket  and  a  typical  sample  mesh,  when 
the  procedure  described  above  has  been  applied  (here,  we  take  6  =  5  and  c2  =  0.002  in  equation 
(29)).  It  should  be  noted  that  with  this  procedure,  one  just  needs  to  generate  a  single  mesh  for 
the  problem  to  be  solved  (the  mesh  of  the  object  with  its  original  shape,  i.e.  the  mesh  of  Figure  3, 
for  example).  One  needs  to  be  careful  that  once  the  points  defining  the  original  object  have  been 
moved  randomly,  the  triangles  sitting  on  the  object  do  not  distort  the  mesh  too  much.  This  can  be 
easily  controlled  by  the  parameter  c  in  equation  (29)  which  is  directly  linked  to  the  amplitude  that 
Sri  can  take.  The  coordinates  of  the  points  which  define  the  finite  element  mesh  only  occur  in  the 
matrices  M,  S,  and  F  of  equation  (16)  and  those  quantities  are  constant  element  by  element.  Since 
those  matrices  appear  as  multiplicative  coefficients  into  Maxwell’s  equations  we  have  transformed 
a  problem  with  an  object  having  a  random  shape  (which  usually  requires  the  generation  of  a  new 
mesh  for  each  sample)  into  an  equivalent  stochastic  problem  having  a  fixed  mesh  but  with  random 
coefficients. 

5  Computation  of  sensitivity  of  the  RCS 

The  (2D)  RCS  under  the  scattering  angle  0,  is  a  function  of  the  following  integral  along  dC  : 


( nx  cos  0  ~h  ny  sin  0)  Ez 

ej2'nf(x  cos  <£+y  sin  <f>)  ^ 


(31) 
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Figure  3:  Original  mesh  for  the  rocket  problem. 


Figure  4:  One  sample  of  a  mesh  for  the  rocket  with  a  random  shape  together  with  its  original  shape. 
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and  normalized  as 


R.CS(<J> )  =  2  Tr/^^i-, 

|E0|2 

where  E0  is  power  of  the  incoming  field.  For  the  2D  calculations,  a  scaling  equal  to  the  wavelength 
A  is  applied  to  the  RCS,  i.e.  the  plots  of  the  RCS  on  a  dB  scale  are  simply 

RCSdB{<}> )  =  10  \ogl0(RCS  ((/>)/ A). 

For  the  RCS  for  3D  scattering  objects,  a  scaling  equal  to  A2  is  applied. 

5.1  Using  Monte-Carlo  simulation 

A  Monte-Carlo  simulation  is  therefore  quite  simple:  one  simply  need  to  generate  (say)  A/  random 
vectors  {£rm}i<m<yv/  as  described  above.  Each  set  of  random  numbers  will  give  a  new  random 
surface  from  which  we  can  compute  M  solutions  of  Maxwell’s  equations.  Those  M  solutions  of 
Maxwell’s  equations,  will  give  M  radar  cross  section  {RCSm} \<m<M  from  which  we  can  compute 
averages  as  follows: 

1  M 

<  RCS  >=jjY  RCS™  (32) 

m—  1 

And  for  the  variance,  we  have 


1  M 

var{RCS )  =  jj  Y  (RCSm)2  ~  <  RCS  >2  .  (33) 

m  =  1 

It  should  be  noted  that  the  advantage  of  the  Monte-Carlo  approach  is  its  simplicity  (it  only  requires 
repetitive  runs  of  an  existing  deterministic  solver),  but  it  is  hard  to  get  accurate  solutions  due 
to  its  slow  rate  of  convergence  0(A/-1/2).  In  the  next  subsection,  we  will  introduce  a  stochastic 
collocation  method  which  has  the  simplicity  of  the  Monte-Carlo  approach  but  with  better  rates  of 
convergence. 

5.2  Computation  of  sensitivity  of  RCS  using  stochastic  collocation 

About  fifty  years  ago.  Stroud  constructed  a  set  of  cubature  points  to  compute  integrals  of  the  form 

/[/l  =  f  f(x)dx.  (34) 

J  l-i.i!n 

This  set  of  cubature  points  based  on  (n  T  1)  points  is  exact  for  polynomials  of  degree  two,  and  the 
approximation  is  written  as 


71+1 


7[/]  -  Y  w«/(Xi)> 

1=1 

where  the  (n  +  1)  cubature  points  x,  =  (x-,x2,  ...,x^)  arc  given  by 

(35) 

(36) 

*  - 
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for  r  =  1, ...,  [n/2],  and  if  n  is  odd,  xj1  =  (-l)^-1)/^.  The  weights  in  (35)  are  all  equal  to  2n/(n+ 1). 
Similarly,  we  have  the  Stroud-3  method  based  on  2 n  points  which  is  exact  for  polynomials  of  degree 
three  : 


2n 

Hf]  -  Ew*/(x<)>  (37) 

i=i 

where  the  2 n  cubature  points  x*  =  (xj,x?,  ...jX-1)  are  now  defined  by 

=  /H(2-1>"),  (38) 

*  - 

for  r  =  l,...,[n/2]  and  if  n  is  odd,  similarly  we  have,  x-1  =  (  —  l)l/\/3.  The  weights  in  (37)  arc  all 
equal  to  2n/2 n.  It  can  be  shown  that  Stroud-2  and  Stroud-3  methods  employ  the  minimal  number  of 
points  for  their  corresponding  integration  accuracy.  Since  Stroud-2  and  Stroud-3  methods  appeared, 
many  other  cubature  formulae  have  been  established  to  compute  various  high-dimensional  integrals. 
In  the  70’s,  Stroud  published  a  book  containing  most  cubature  formulae  known  at  that  time.  This 
extensive  work  wets  then  continued  by  Cools  in  a  series  of  two  papers.  The  idea  of  the  stochastic 
collocation  method  is  based  on  polynomial  interpolations  in  the  multi-dimensional  random  space. 
We  assume  that  Lagrange  polynomials  based  on  Stroud's  cubature  points  exist,  and  we  express  the 
RCS  using  the  Lagrange  interpolation  polynomials,  which  gives  (for  Stroud-2  cubature  points) 


n+1 

RCS(x)  =  RCSiLi(x),  (39) 

1=1 

where  {Li}i<i<n+i  are  (n  +  l)-variate  Lagrange  polynomials  based  on  points  {xi}i<j<n+i  of  the 
cubature  formula  (35).  We  note  that  by  construction  of  Lagrange  polynomials,  we  have  Lt(xj)  =  6tJ 
and  therefore  RCSi  =  RCS(xt).  The  radar  cross  sections  RCS(xt)  can  be  easily  computed  as 
follows: 

•  for  each  cubature  point  x?  generate  a  new  surface  with  <Srt  =  PD1/2G,  and  G,  =  y/3 x* 

{ i  =  1, n  T  1) 

•  compute  the  corresponding  solution  of  Maxwell’s  equations 

•  compute  the  corresponding  radar  cross  section  RCS(xj) 

The  expression  of  the  RCS  is  now  available  under  the  form  (39)  and  we  will  show  that  statistical 
quantities,  e.g.  mean  and  variance,  can  easily  be  computed.  By  taking  the  average  of  equation  (39) 
and  evaluating  the  multi-dimensional  integral  with  Stroud-2  cubature  formula  (35),  we  get 

1  ri+ 1  r 

<  RCS(x)  >=  —  y  RCSi  /  Lt(x)dx 

2n  fr{  J\-i,i]* 

I  n  +  1  n- 1-1 

^  i=i  j=i 
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n+1  n+ 1 

= 

Z=1  j=l 

1  n+l 

Z  Z=1 

Similarly,  for  the  variance,  we  have 

i  n+l 

var{RCS(x))  =  —  Y,  UiRCS?-  <  RCS(x)  >2  .  (40) 

2  <=i 

The  same  procedure  can  be  used  with  Stroud-3  cubature,  and  in  that  case,  2 n  realizations  of  the 
RCS  corresponding  to  the  2 n  cubature  points  of  equation  (37)  will  have  to  be  computed.  When 
the  random  space  is  one-dimensional  this  procedure  can  also  be  repeated  for  the  Legend  re- Gauss- 
Lobatto  quadrature. 

Alternatively,  we  could  have  used  orthogonal  polynomials  to  express  the  RCS.  i.e.  equation  (39) 
would  have  to  be  replaced  by 


n+l 

RCS(x)  =  Y,  RCSi$i(x),  (41) 

1  =  1 

where  4>i(x)  are  orthogonal  polynomial  on  [— l,l]n.  The  expression  (41)  is  usually  referred  as 
” polynomial  chaos  expansion”  and  was  used  in  a  number  of  mechanical  problems.  The  problem  of 
this  approach  lies  in  the  difficulty  to  compute  in  an  efficient  way  the  coefficients  RCSi. 

5.3  Reduced  random  space 

For  the  stochastic  collocation  method,  we  have  seen  that  there  is  a  close  relation  between  the  size  of 
the  random  space  n  and  the  number  of  calls  to  the  deterministic  Maxwell  solver  (i.e.  we  need  (71  + 1) 
calls  for  Stroud-2  and  2 n  calls  for  Stroud-3).  A  way  of  reducing  the  CPU  cost  of  the  method  would 
be  to  reduce  the  size  of  the  random  space  n.  This  is  possible,  up  to  a  certain  extent,  depending  on 
the  covariance  matrix  K  defined  by  equation  (29).  Let  us  consider  the  two  extreme  cases: 

First,  we  assume  that  6  — >  0  and  c  =  1  in  equation  (29);  in  this  case,  K  =  I  is  the  identity 
matrix  and  P  =  D  =  I.  In  other  words  we  will  need  n  independent  and  identically  distributed 
random  variables  to  describe  the  process  since  all  points  of  the  object  are  uncorrelated. 

Now  we  assume  that  b  — ►  00  and  c  =  1  in  equation  (29);  in  this  case,  all  variables  arc  fully 
correlated.  We  have  K{j  =  Pij  =  1  for  i,j  =  1,  ...,n;  D\\  =  1  and  Da  =  0  for  i  =  2,  ...,n.  Thus,  all 
n  points  of  the  sample  5r  ={5ri,Sr21  ...,(5rn)7  =  PD1/2G  are  moved  with  the  same  amplitude,  i.e. 
Sr\  =  6r2  =  ...  =  6rn.  This  means  that  n  —  1  would  suffice  to  describe  the  problem.  In  practice,  we 
are  between  those  two  extreme  cases  and  the  size  of  the  random  space  can  be  given  by  the  number 
of  significant  eigenvalues  of  D  denoted  by  m  (with  rn  <  n).  Then,  the  modified  algorithm  proceeds 
as  follows: 

•  Compute  the  n  eigenvalues  {A*}i<j<n  of  the  matrix  K  and  select  the  m  most  significant  ones 
satisfying  >  TOL ,  where  TOL  is  some  small  number.  The  integer  rn  will  be  the  size  of 

the  random  space  needed  to  move  the  points  of  the  objects  randomly. 
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•  Wc  denote  by  P  and  D  the  matrices  of  size  n  x  m  constructed  from  the  restriction  of  the 
matrices  P  and  D  where  we  have  only  kept  the  most  significant  eigenvalues  and  their  associated 
eigenvectors. 

•  Since  the  dimension  of  the  random  space  is  reduced  to  m  <  n,  the  collocation  points  x 
are  m- dimensional  vectors  based  on  quadrature  points  for  the  computation  of  7n-dimensional 
integrals.  For  Stroud-2,  we  have  (m  +  1)  points  and  for  Stroud-3,  it  is  2m. 

•  The  vector  Sr  used  to  move  the  mesh  is  now  given  by  Sr  =(5ri,  <5r2, £rn)T  =  PD]/2G, 
where  G  =  \/3x. 


6  A  few  examples 

We  now  consider  a  square  of  width  0.8A  whose  mesh  consists  of  840  elements.  Again,  A  is  the 
wavelength  of  the  incident  plane  wave.  The  dimension  of  the  full  random  space  is  44,  which  means 
that  the  square  is  formed  by  44  points  of  the  finite  clement  mesh.  The  parameters  6  and  c2  in 
equation  (29)  are  taken  equal  to  5  and  0.001,  respectively.  This  choice  of  c  corresponds  to  the 
radial  randomness  having  a  standard  deviation  of  about  0.032 A.  By  choosing  6  =  5,  the  covariance 
matrix  K  defined  by  (29)  is  more  diagonal  dominant  than  it  was  in  the  previous  example,  and  the 
number  of  most  significant  eigenvalues  is  reduced  to  m  =  16.  Due  to  the  geometrical  singularities 
in  the  corners  of  the  square,  this  is  a  much  harder  problem  to  solve  than  the  cylinder  problem. 
This  is  illustrated  on  Figure  5  and  6,  where  we  show  the  mean  and  the  variance  of  the  RCS  for  the 
Stroud-2,  the  Stroud-3  and  the  Monte-Carlo  methods  using  250  samples.  We  can  see  on  Figure  5 
that  the  mean  converges  to  the  same  values  for  the  three  methods.  However,  for  the  variance,  the 
Stroud-3  and  the  Monte-Carlo  methods  give  similar  results  but  the  Stroud-2  method  gives  results 
which  have  not  converged,  especially  in  the  sidebands.  This  is  because  Stroud-2  can  only  integrate 
exactly  multi-variate  polynomials  of  degree  two  at  most,  and  the  exact  solution  of  this  example 
cannot  be  accurately  represented  by  such  polynomials.  On  the  other  hand,  Stroud-3  which  can 
integrate  exactly  multi-variate  polynomials  of  degree  three  at  most  does  a  better  job.  In  order  to 
emphasize  the  cost  saving  of  Stroud-3  method  over  Monte-Carlo  method,  we  have  represented  the 
variance  of  the  RCS  for  Stroud-3,  Monte-Carlo  with  32  samples  (which  has  the  same  CPU-cost  as 
the  Stroud-3  method)  and  Monte-Carlo  with  250  samples  in  Figure  7.  Wc  see  that  as  wc  increase 
the  number  of  samples,  the  Monte-Carlo  solution  converges  slowly  to  the  Stroud- 3  solution.  Figure 
8  shows  the  average  of  the  RCS  and  the  possible  variations  around  its  average  value  obtained  with 
Stroud-3  cubature.  We  can  see  that  like  for  the  cylinder  problem,  the  uncertainty  in  shape  affects 
the  RCS  mostly  in  the  sidebands. 

6.1  General  shape 

As  a  last  shape,  wc  consider  the  (simple)  rocket  shown  in  Figure  3  which  is  A  long  and  the  body  is 
0.25A  wide.  The  mesh  consists  of  1465  elements  and  the  dimension  of  the  full  random  space  is  45 
(i.e.  the  rocket  is  formed  by  45  points  of  the  finite  element,  mesh).  As  for  the  square  problem,  the 
parameters  6  and  c2  in  equation  (29)  are  taken  equal  to  5  and  0.001,  respectively  and  the  size  of 
the  random  space  can  be  reduced  to  m  =  14.  This  choice  of  c  corresponds  to  the  radial  randomness 
having  a  standard  deviation  of  about  0.032A.  Similarly  to  the  square  problem,  good  convergence 
of  the  average  RCS  is  obtained  both  for  Stroud-2,  Stroud-3  and  Monte-Carlo  simulations.  However 
for  the  variance  of  the  RCS,  only  Stroud-3  and  Monte-Carlo  with  enough  samples  (i.e.  more  than 
250)  give  good  results.  Figure  9  shows  the  average  of  the  RCS  and  the  possible  variations  around 
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Figure  5:  Comparison  of  the  average  of  the  RCS  for  the  square  problem  for  Stroud-2,  Stroud-3  and 
Monte-Carlo  methods. 


Figure  6:  Comparison  of  the  variance  of  the  RCS  for  the  square  problem  using  Stroud-2,  Stroud-3 
and  Monte-Carlo  methods. 
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Figure  7:  Comparison  of  the  variance  of  the  RCS  for  the  square  problem  using  Stroud-3  method 
and  Monte-Carlo  method  with  32  and  250  samples. 


Figure  8:  RCS  for  the  square  problem.  Results  are  shown  with  the  mean  RCS  as  well  as  ±  one 
standard  deviation. 


17 


Figure  9:  RCS  for  the  rocket  problem.  Results  are  shown  with  the  mean  RCS  as  well  as  ±  one 
standard  deviation. 

its  average  value  obtained  with  Stroud-3  cubature.  For  this  case,  the  uncertainty  in  shape  affects 
the  RCS  both  in  the  middle  part  and  in  the  sidebands. 

6.2  Higher  frequency  cases 

As  a  last  numerical  experiment  in  two  spatial  dimensions,  we  increase  the  frequency  /  of  the  source 
term  (??)  up  to  /  =  3.  This  causes  the  radius  of  the  cylinder  to  increase  to  3 A  and  the  size  of  the 
rocket  to  increase  to  3 A  long  and  0.75 A  wide  in  the  body.  For  these  cases  the  radial  randomness  has 
a  standard  deviation  of  about  0.095A.  For  all  the  previous  numerical  experiments,  we  have  used 
4th  order  elements  in  each  triangle  of  the  finite  element  grid.  Due  to  the  higher  frequency,  all  the 
results  shown  in  this  section  are  obtained  with  6th  order  elements  to  get  converged  results  in  the 
physical  space  (i.e.  the  value  of  p  defined  below  equation  (11)  is  increased  to  p  =  6).  Figures  10 
and  1 1  show  the  average  of  the  RCS  and  the  possible  variations  around  its  average  value  obtained 
with  Stroud-3  cubature  for  the  cylinder  and  the  rocket  problem,  respectively.  The  form  of  the 
uncertainty  remains  the  same  as  in  the  previous  numerical  examples,  the  only  difference  being  the 
frequency  of  the  source  term.  By  comparing  Figure  ??  with  Figure  10  and  Figure  9  with  Figure 
11,  we  see  that  as  we  increase  the  frequency,  the  regions  where  the  RCS  has  the  greatest  variance 
remains  more  or  less  the  same.  In  this  numerical  experiment,  convergence  of  the  average  of  the  RCS 
and  the  variance  of  the  RCS  is  good  for  the  cylinder  problem,  for  both  Stroud  and  Monte-Carlo 
methods.  However,  for  the  rocket  problem  only  Stroud-3  and  Monte-Carlo  with  enough  samples 
give  converged  results  for  the  variance. 

6.3  Sphere 

For  a  first  example  with  three  spatial  dimensions  we  consider  the  scattering  of  a  transverse  magnetic 
plane  wave  from  a  PEC  sphere.  We  assume  the  sphere  has  a  uniform  random  radius  in  the  interval 
[0.9A,  1.1  A],  where  A  is  the  wavelength  of  the  incident  field.  The  use  of  one  random  variable  is 
not  a  limitation  of  the  method  but  is  done  for  logistical  reasons  to  have  reasonable  computation 
times.  Since  this  experiment  has  only  one  random  dimension,  a  sixth  order  Legendrc-Gauss-Lobatto 
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Figure  10:  RCS  for  the  cylinder  problem  at  higher  frequency.  Results  are  shown  with  the  mean 
RCS  as  well  as  ±  one  standard  deviation. 


Figure  11:  RCS  for  the  rocket  problem  at  higher  frequency.  Results  are  shown  with  the  mean  RCS 
as  well  as  ±  one  standard  deviation. 
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Figure  12:  One  sample  of  a  surface  mesh  for  the  sphere  with  a  random  radius. 


Figure  13:  RCS  for  the  sphere  problem.  Results  are  shown  with  the  mean  RCS  as  well  as  ±  one 
standard  deviation. 

quadrature  is  used  for  collocation  in  the  random  space.  For  the  spatial  discretization  we  use  fourth 
order  elements  and  a  sample  mesh  is  presented  in  Figure  12  which  is  restricted  to  the  surface  of 
the  sphere.  Figure  13  shows  the  average  of  the  RCS  and  the  possible  variations  around  its  average 
value.  The  uncertainty  in  the  radius  of  the  sphere  affects  the  RCS  mainly  in  the  sideband. 

6.4  Three-dimensional  rocket 

For  the  second  experiment  with  three  spatial  dimensions  we  consider  the  scattering  of  a  transverse 
magnetic  plane  wave  from  a  PEC  rocket.  The  orientation  of  the  scattering  plane  and  the  scattering 
angle  9  is  is  given  in  Figure  14  along  with  a  geometric  description  of  the  rocket.  The  scattering  angle 
of  the  incident  field  is  assumed  to  be  random  with  uniform  distribution  in  the  interval  [10°,  20°].  As 
in  the  case  for  the  sphere,  the  use  of  one  random  variable  is  not  a  limitation  of  the  method  but  is 
done  for  logistical  reasons  to  have  reasonable  computation  times.  A  fourth  Legendre-Gauss-Lobatto 
is  used  for  collocation  in  the  one-dimensional  random  space.  For  this  calculation  the  physical  space 
is  discretized  with  degree  five  polynomials  in  each  clement.  Figure  15  shows  the  mesh  used  restricted 
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Figure  14:  The  orientation  of  the  rocket  with  respect  to  the  scattering  plane  is  given  in  (a).  Here  0 
is  the  scattering  angle.  The  top  view  presents  a  projection  of  the  rocket  onto  the  scattering  plane. 
The  geometric  description  of  the  three-dimensional  rocket  is  given  in  (b).  Here  A  is  the  wavelength 
of  the  incident  field. 
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Figure  15:  The  surface  mesh  for  the  three-dimensional  rocket. 

to  the  surface  of  the  rocket  and  Figure  16  shows  the  average  of  the  RCS  and  the  possible  variations 
around  its  average  value.  The  uncertainty  in  the  direction  of  the  incident  field  affects  the  RCS 
mainly  near  the  local  maxima  and  minima  points  of  the  RCS. 
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Figure  16:  RCS  for  the  three-dimension  rocket  problem.  Results  are  shown  with  the  mean  RCS  as 
well  as  ±  one  standard  deviation. 
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10  Transition 

As  outlined  in  the  original  proposal,  it  was  planned  that  some  of  the  techniques  being  developed  in 
this  effort  be  implemented  into  the  TEMPUS  CEM  code,  developed  and  maintained  by  HyperCornp, 
Inc  and  used  extensively  for  RCS  prediction  and  modeling  by  DoD  contractors  and  laboratories. 

This  has  been  accomplised  and  we  stay  in  continued  close  contact,  The  PI  has  visited  HyperCornp 
several  times  during  the  period  to  report  on  progress  and  discuss  ways  of  implementing  the  developed 
techniques  into  TEMPUS  which  now  has  elementary  support  for  uncertainty  quantification  During 
the  effort,  we  have  also  stayed  in  close  contact  with  WPWFB  (Dr  Camberos)  to  update  him  on  status 
and  enable  a  possible  future  transition  of  the  technology  to  their  internal  simulation  framework. 
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